﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Nature;
using Nature.Numerics;
using System.Diagnostics;
using Nature.Numerics.Basic;

namespace Nature.Toolkit.GEqu
{
    public class TPIdealGasChemicalEquilibriumState : IChemicalEquilibriumState, ISoneFunction, ISoneApproximation
    {
        #region Private Fields

        private SoneNewton<TPIdealGasChemicalEquilibriumState> m_newton;
        private bool m_equilibrated = false;
        private double? m_t;
        private double? m_p;
        private bool m_compositionInitialized = false;

        /// <summary>
        /// Species mole fraction / variable (Lagrangian coefficient) derivatives matrix
        /// </summary>
        private double[][] m_DxDa;
        /// <summary>
        /// Element switches
        /// </summary>
        bool[] m_elementSwitches;
        /// <summary>
        /// Species mole fractions evaluated based on the current variable (Lagrangian coefficients)
        /// </summary>
        private double[] m_x;
        /// <summary>
        /// Species mole fractions evaluated based on the selected variable (Lagrangian coefficients).        
        /// </summary>
        /// <remarks>
        /// Linear approximation base point
        /// </remarks>
        private double[] m_x0;
        /// <summary>
        /// Variable increment (Lagrangian coefficients)
        /// </summary>        
        private double[] m_da;
        /// <summary>
        /// Linear approximation base point (Lagrangian coefficients).   
        /// </summary>
        private double[] m_a0;
        /// <summary>
        /// Logarithm of the gas pressure (problem parameter) reduced with the atmospheric pressure
        /// </summary>
        private double m_logP;
        /// <summary>
        /// Reference elements molar composition [mol/kg] (problem parameter)
        /// </summary>
        private double[] m_n0;

        /// <summary>
        /// Elements molar composition [mol/kg]
        /// </summary>
        private double[] m_n;

        private M1KgUtil m_m1kgUtil;

        private double[] m_variable;

        private double[] m_gibbs;

        #endregion

        #region Private Methods

        private void Equilibrate()
        {
            if (m_equilibrated) { return; }
            if (false == m_t.HasValue) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullTemperatureException(); }
            if (false == m_p.HasValue) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullPressureException(); }
            if (false == m_compositionInitialized) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullCompositionException(); }
            object anchor = null;
            m_newton.Solve(this.m_variable, ref anchor);
            this.Exact_x_using_a(m_x, m_variable);
            Debug.Assert(Math.Abs(m_x.Sum() - 1.0) < 1.0e-4, "|sum(x_k) - 1| < 1.0e-4");
            m_equilibrated = true;
        }

        #endregion

        #region Constructors

        public TPIdealGasChemicalEquilibriumState(ModelObject model)
        {
            this.Model = model;

            m_newton = new SoneNewton<TPIdealGasChemicalEquilibriumState>(this)
            {
                MaxNumberOfLinearApproximationReuse = 1
            };
            m_newton.SetAbsoluteConvergenceRadius(10);
            m_newton.SetRelativeConvergenceRadius(1.0e-1);
            m_newton.SetAbsoluteTolerance(1.0e-5);
            m_newton.SetRelativeTolerance(1.0e-5);
            m_newton.SetUpperLimit(+1000);
            m_newton.SetLowerLimit(-1000);

            m_gibbs = new double[model.Nsp];
            m_variable = new double[model.Nel];

            m_m1kgUtil = new M1KgUtil(model);
            m_DxDa = Numerics.Utilities.CreateMatrix<double>(model.Nsp, model.Nel);
            var elspNormalized = model.ElspNormalizedMatrix;

            m_elementSwitches = new bool[Nel];
            for (int iel = 0; iel < Nel; ++iel)
            {
                m_elementSwitches[iel] = true;
                var curr = elspNormalized[iel];
                for (int l = iel + 1; l < Nel; ++l)
                {
                    var next = elspNormalized[l];
                    if (curr.IsSame(next))
                    {
                        m_elementSwitches[iel] = false;
                    }
                }
            }

            m_a0 = new double[model.Nel];
            m_da = new double[model.Nel];
            m_n0 = new double[model.Nel];
            m_n = new double[model.Nel];

            m_x0 = new double[model.Nsp];
            m_x = new double[model.Nsp];
        }
        #endregion

        #region Public Properties

        public ModelObject Model { get; private set; }

        /// <summary>
        /// Number of elements
        /// </summary>
        public int Nel { get { return Model.Nel; } }

        /// <summary>
        /// Number of species
        /// </summary>
        public int Nsp { get { return Model.Nsp; } }
        #endregion

        #region Private Methods

        /// <summary>
        /// Computes mole fractions vector for the given variables vector based on the exact operator form
        /// </summary>
        /// <param name="x">mole fractions</param>
        /// <param name="a">variable</param>
        private void Exact_x_using_a(double[] x, double[] a)
        {
            var spel = this.Model.SpelMatrix;
            double pmin = Math.Log(1.0e-30);
            double pmax = Math.Log(1.0e+3);
            double p;

            for (int isp = 0; isp < Nsp; ++isp)
            {
                p = m_gibbs[isp] + m_logP;
                for (int iel = 0; iel < Nel; ++iel)
                {
                    p += spel[isp][iel] * a[iel];
                }
                p = (-p);
                if (p < pmin) p = pmin;
                else if (p > pmax) p = pmax;
                x[isp] = Math.Exp(p);
            }
            //Debugger.Break();
        }

        /// <summary>
        /// Computes mole fractions vector for the given variables vector based on the approximate operator form
        /// </summary>
        /// <param name="x">mole fractions</param>
        /// <param name="a">variable</param>
        private void Approximate_x_using_a(double[] x, double[] a)
        {
            for (int iel = 0; iel < Nel; ++iel)
            {
                m_da[iel] = (a[iel] - m_a0[iel]);
            }
            for (int isp = 0; isp < Nsp; ++isp)
            {
                double dx = 0.0;
                for (int iel = 0; iel < Nel; ++iel)
                {
                    dx += m_DxDa[isp][iel] * m_da[iel];
                }
                x[isp] = m_x0[isp] + dx;
            }
        }
        #endregion

        #region IFunction Members

        /// <summary>
        /// Number of variables
        /// </summary>        
        [DebuggerNonUserCode()]
        int ISoneFunction.Nva { get { return this.Nel; } }

        /// <summary>
        /// Computes the function vector for a given variable vector using the exact operator form
        /// </summary>
        /// <param name="a">variable vector</param>
        /// <param name="eq">function vector</param>
        public void Compute(double[] a, double[] eq)
        {
            Exact_x_using_a(m_x, a);
            ComputeUsingX(m_x, eq);
            for (int iel = 0; iel < Nel; ++iel)
            {
                if (m_elementSwitches[iel] == false)
                {
                    eq[iel] = a[iel];
                }
            }

        }

        /// <summary>
        /// Computes the function vector based on the given mole fraction values
        /// </summary>
        /// <param name="x">mole fractions</param>
        /// <param name="eq">function vector</param>
        private void ComputeUsingX(double[] x, double[] eq)
        {
            double xsum = 0.0;
            for (int isp = 0; isp < Nsp; ++isp) { xsum += x[isp]; }
            m_m1kgUtil.ComputeElementMolesPerMass(m_n, x);
            for (int iel = 0; iel < Nel; ++iel)
            {
                eq[iel] = m_n[iel] * xsum - m_n0[iel];
            }
        }

        #endregion

        #region ISoneApproximation Members

        public void ConfigureApproximation(double[] a0)
        {
            var spel = Model.SpelMatrix;
            m_a0.Assign(a0);
            Exact_x_using_a(m_x0, m_a0);
            for (int isp = 0; isp < Nsp; ++isp)
            {
                for (int iel = 0; iel < Nel; ++iel)
                {
                    double akl = spel[isp][iel];
                    double x = Math.Max(1.0e-8, m_x0[isp]);
                    m_DxDa[isp][iel] = (x * (-akl));
                }
            }
        }


        int ISoneApproximation.Nva
        {
            get { return Model.Nel; }
        }


        void ISoneApproximation.Approximate(double[] a, double[] f)
        {
            Approximate_x_using_a(m_x, a);
            this.ComputeUsingX(m_x, f);
            for (int iel = 0; iel < Nel; ++iel)
            {
                if (m_elementSwitches[iel] == true) { continue; }
                f[iel] = a[iel];
            }
        }

        #endregion

        #region IChEqCalculator Members

        double IChemicalEquilibriumState.T
        {
            get
            {
                if (m_t.HasValue) { return m_t.Value; }
                throw ChemicalEquilibriumStateNotInitializedException.CreateNullTemperatureException();
            }
            set
            {
                if (m_t != value)
                {
                    m_t = value;
                    m_equilibrated = false;
                    this.Model.IdealGasSpeciesThermoPotentials.G(m_t.Value, m_gibbs);
                }
            }
        }

        double IChemicalEquilibriumState.P
        {
            get
            {
                if (m_p.HasValue) { return m_p.Value; }
                throw ChemicalEquilibriumStateNotInitializedException.CreateNullPressureException();
            }
            set
            {
                if (m_p != value)
                {
                    m_p = value;
                    m_equilibrated = false;
                    m_logP = Math.Log(m_p.Value / Constants.Atmosphere);
                }
            }
        }

        ReadOnlyArray<double> IChemicalEquilibriumState.X
        {
            get
            {
                Equilibrate();
                return m_x;
            }
            set
            {
                m_compositionInitialized = true;
                if (value.Length != Nsp) { throw new ArgumentException("value.Length != Nsp"); }
                Debug.Assert(value.All(i => i >= 0.0 && i <= 1.0));
                Debug.Assert(Math.Abs(1.0 - value.Sum()) < 1.0e-8);
                bool changing = false;
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    changing |= (m_x[isp] != value[isp]);
                    m_x[isp] = value[isp];
                }
                if (changing)
                {
                    m_equilibrated = false;
                    //m_m1kgUtil.ComputeElementMolesPerMass(m_n, value);
                    this.Model.CalcElementAmountPerMass(m_n0, value);
                }
            }
        }

        public void DropState()
        {
            m_variable.Assign(0.0);
            m_equilibrated = false;
            m_p = null;
            m_t = null;
            m_compositionInitialized = false;
        }

        #endregion
    }
}